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The equations of motion of a secularly precessing elliptical 
orbit 
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ABSTRACT 

The equations of motion of a secularly precessing ellipse are developed using time as 
the independent variable. The equations are useful when integrating numerically the 
perturbations about a reference trajectory which is subject to secular perturbations 
in the node, the argument of pericenter and the mean motion. Usually this is done 
in connection with Encke's method to ensure minimal rectification frequency. Similar 
equations are already available in the liter ature, but they are either given based on the 
true anomaly as the independent variable ()Kvner fc Bennetdll966h . or in mixed mode 
with respect to the time through the use of a supporting equation to track the anomaly 
(|Escobal| [l965') . The equations developed here form a complete and independent set of 
six equations in the time. Reformulations both of Escobal's and Kyner and Bennett's 
equations are also provided which lead to a more concise form. 

Key words: method: analytical - celestial mechanics 



1 INTRODUCTION 

In many cases one is confronted with the evaluation of the spectral characteristics of specific perturbation accelerations on 
the orbital motion of planetary satellites. In several instances analytical evaluations can be provided. However, this is not 
always possible, and occasionally one may desire verification of the analytical results. Since the analytical development of 
perturbations is always based on some reference orbit, it becomes of interest in numerical verification work to be able to 
reproduce the same reference orbit. The simplest such orbit is an invariable ellipse, and Encke's method can then be adopted 
for the numerical integration of perturbations about this nominal orbit. Encke's method can be generalized to include non- 
Keplerian reference trajectories. The method can then be formulated in terms of the differential equation of the reference orbit 
r* (t) and the equation of the perturbation of the relative position Sr = r (t) — r* (t) of the trajectory r (t) . These equations 
are 



d 2 r* 

~dlF 

d 2 8r 
dt 2 



—^r + s (r 



r*,t), 

{f(q)r* + [l + f (q)] Sr}+g(r,r,r*,r*,t), 



(1) 
(2) 



where ii = G (mi + iri2) is the gravitational parameter of two masses mi and m,2, s(r,r,t) is the disturbing acceleration 
which generates the non-Keplerian nominal trajectory, while g is the complementary disturbing acceleration. This is defined 
as the difference between the actual disturbing acceleration p (r, r, t) acting on the body whose orbit is being propagated and 
the disturbing acceleration of the nominal trajectory, 



g(r,r,r*,r*,t) =p(r,r,t) - s(r*,r*,t). 

Encke's parameter q is defined |Battinll 19871 ') as 

(2r* + Sr) ■ Sr 
(r* + Sr) ■ (r* + 5r) ' 

and the auxiliary function / (q) has the form 



(3) 



(4) 
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3 1 + g + g ,^ 
l + (l + q) 3/2 

Note that in the unperturbed case (s — p = 0) Encke's equations |T} and © are exactly equivalent to the Two-Body equations 
of motion of the two trajectories. Of course, in this case, their use is limited to non-osculating initial conditions. In general, 
however, the initial conditions for the perturbed case are such that Sr (to) = 5r (to) = 0. Also note that Encke's equations are 
completely equivalent to the original perturbation problem stated for trajectories r* (t) and r (t) . In analytical developments 
the approximation is usually made, however, of evaluating the disturbing acceleration p on the nominal trajectory, that is, of 
repl acing q (r,r,r* ,r* , t) wi th g (r* , r* ,t) = p (r* , r* ,t) — s (r* , r* ,t) . 

iKvner fe Bennett] (1 19661 ) (hereafter K&B) proposed a modification to Encke's method which replaces the invariable 
ellipse with a precessing ellipse that includes the first-order secular effects of the Earth's oblateness. This has the desirable 
effect of limiting the frequency of rectifications necessary to contain within specified bounds the unavoidable divergence of 
the perturbed orbit from the reference orbit. The K&B formulation is given using the true anomaly as the independent 
variable . For compa rison with analytical theories based on time as the independent varia ble, like Kaula' s linear satellite 
theory jKaulal 19661). the K&B approach needs to be ref ormulated. This has be en done by lEscobal (|l965T ). More recently 
iLundberg et all l|200d ). on the basis of previous work by ILundberg et all |l99ll ). have introduced an extension of Encke's 
method for application to long-arc orbit determination which uses a precessing and librating ellipse of variable shape based on 
the Escobal model. The behaviour of the K&B and the Escobal reference orbits are similar, but the K&B approach contains 
some periodic terms in addition to the purely secular terms of the time-wise approach of Escobal. Comparison with Kaula's 
theory requires a time-wise approach, but implementation of Escobal's original equations is cumbersome since they require 
the use of a supporting equation for the true anomaly, also a feature of the K&B approach. It is then desirable to investigate 
the possibility of improving on the available representations of the secularly precessing ellipse. In the following we will review 
these classical methods and provide for each an alternative formulation. Finally we will develop a novel formulation with 
the time as independent variable, but distinctly different from Escobal's, which leads to a vector differential equation of the 
second order, whose coefficients are functions of the dynamical state variables only. 



2 THE SECULARLY PRECESSING ELLIPSE 

The secularly precessing ellipse (SPE) can be defined through the secular rates of the angular elements either as a function of 
the true anomaly or as a function of the mean anomaly. Here we define the SPE for both cases providing analytical expressions 
for the secular rates of the angular elements due to the first zonal harmonic of the gravity field of the central body. 

2.1 True anomaly / as independent variable 

Given the gravitational parameter fxo = G (mi +7712) , let the orbital elements at epoch be (ao, eo, io, Ho, ^0, Mo) where the 
associated Keplerian mean motion is no = f/io/«o) 2 ■ in the presence of a second degree zonal harmonic potential J2 = — C20, 
the secular variations of the elements can be given as a function of the true anomay / as 

a = ao, fl(f) = r(/-/o) + n , 

e = e , w (/) = T) (/ - f ) + wo, (6) 

i = i , M (t) = n(t-t ) + Mo, 

where the perturbed secular mean motion is given by 

n = n (1 - 7) , (7) 
and the constant secular rates r and n due to the first zonal harmonic are given bv lSterne] l|l960l l IKvner fc Bennett] (|l966h 

dO. 3 J 2 ^ 2 



df ~ 2(l- e 2) 



duo 3 J2 

V = 



df 4 (1 



2 (^) cosi > ( 8 ) 
(^) 2 (4 -5 sin 2 *), (9) 



with 

7 = "^ j2 (?) 2 (^) 3 [ 1_3sin2isin2(aJ(, + /o) ]- (10) 
where ro = a (l — e 2 ) / (1 + ecos/o) . Note that the phoronomic elements a, e and i do not show secular behavior. 

2.2 Mean anomaly M as independent variable 

If we choose time t as the independent variable, then the secular rates of the orbital elements assume the expressions 
a = ao, f2 (t) = Cl (t — to) + Oo, 
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e = e , w (*) = w(< — to) + wo, (11) 
i = i 0) M (t) = n (i - to) + M , 

where the constant secular rates D. and u due to the first zonal harmonic are given bv lKaulal (|l96rj ) 

a 3 n C 2 o /a e \ 2 . 

S2 = ~ — cosj, (12) 



3 noC2o 



^) 2 (l_5cos 2 i), (13) 



4(1- 

and the perturbed mean motion n is 

n = 710(1-7), (14) 
with 

3 C20 



7 ; 



4(1 



e2) 3/2 



(^) (l-3cos 2 i). (15) 



2.3 Mapping the SPE to cartesian state 

Equations © and equations pip provide at any time t the secular evolution of the orbit elements of the reference orbit. As 
noted above, the two reference orbits are not exactly the same since the true anomaly is a periodic function of time. Given 
the initial Keplerian state vector (ao, eo, io, flo , 0J0, Mo) at epoch to, and the secular rates (J8]l and (|9} and the constant (|10f) 
in the case of true anomaly as independent variable, or the secular rates ()12[1 and (|13p and the constant (|15|l if the mean 
anomaly is adopted as the independent variable, the nominal or reference trajectory r*(t) — hereafter simply indicated with 
r(t) for conciseness — is completely specified at any time t. The corresponding Cartesian state r(t),v(t) can be computed by 
the following procedure 

(i) Compute the nominal true anomaly / by solving first the modified Kepler's equation 

n(t- t ) + Mo = E - e sin E, (16) 

for the eccentric anomaly E and then by transforming to / via the usual Gauss formula 



-W^S-!- (I7) 

The modified mean motion h is either n from equation {7) or n from equation (|14[1 . depending respectively on whether the 
true anomaly or the time is assumed as the independent variable. 

(ii) Compute the radius vector r and the radial velocity v r by means of 



a(l-e 2 ) haesinf 
1 + e cos / ' r VI ~ e 2 ' 

(iii) Finally position r(t) and velocity v(t) are given by 



(18) 



r(t) = R(0,j> + f)r r tn, (19) 

v(t) = R(n,i,uj + f)r rtn + R(n,i,ij + f)r rtn , (20) 

where 

rrtn = I I , r rtn =10), r rtn = I I , (21) 






the vector r rtn having been defined for later use, and the rotation matrix R (fi, i, u + /) from the orbital reference frame (the 
RTN frame, defined by the local radial, transverse, and normal directions) to the inertial reference frame is given by 

R («, i, u) + f) = D (fi) C (i) B (a; + /) , (22) 

with the elementary rotations defined as 

(cos Q — sin Q \ 
sinfi cosfi , (23) 
1/ 

/ 1 \ 

C (i) = I cos i — sin i , (24) 
V sini cosi / 
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/ cos(w + /) -sin(w + /) \ 
B (u + /) = sin (w + /) cos (u + f) , (25) 
V 1 / 

and where the node f2 and the argument of pericenter ui have been updated using equations (JS| or fTTJ. 



2.4 Comparison of the reference orbits 

The two reference trajectories defined above are clearly different. The main difference between the two formulations is due 
to the equation of the center, that is, to the periodic nature of the difference between the true and the mean anomaly. 
Fundamentally this implies that the SPE in true anomaly oscillates periodically with respect to the SPE in mean anomaly. 



3 THE SPE IN THE TRUE ANOMALY 

This Section provides the formulation of the secular ly precessing ellipse with the true anomaly as the independent variable. 
We first review the equations of motion developed bv lKvner fc Bennett! (|l966T l and then provide a more concise reformulation. 



3.1 The Kyner and Bennett formulation 

The equations of motion can be derived easily by differentiating twice equation (|19[) with respect to t and taking into account 
equation ©. We omit the derivation and give only the equations of motion which describe the motion on SPE. The formulation 
is essentially the same as the original formulation of Kyner and Bennett, that is 

.. u (l-7) 2 , 2 (l + ecos/) f 2 d 2 D„„ „ ,„ N dD„dB , 2 „ \ „1 ,„„x 

r + V = "° (1 - ,a 11 [r V CB + 2T{1 + "> *i C *T - +2 ^ R \ ^ (26) 

the only differences being that we have factored (f — 7) 2 on the right hand side and summed the Newtonian term with the 
last term on the right hand side of the original K&B equation (12). Note that here we have introduced the argument of 
latitude u — uj + f. The true anomaly appears explicitly in the equations of motion ()26|) through the rotation matrices B and 
D (equations [25] and [23] respectively). It is thus necessary to compute / either from the analytic procedure outlined previously, 
or by numerically integrating the auxiliary equation 

f= (1 _ e2) a/2 ( 1 + ecos /) • ( 27 ) 



simultaneously with (|26l) . Note that although the perturbed mean motion n appears, this equation is independent of the 
radius vector r. The term 1 + e cos / in equations p6[) cannot be replaced with its Keplerian equivalent a(l — e 2 )/r unless the 
radius vector appearing at the denominator be computed by its Keplerian definition — which means going back to using the 
true anomaly again, or using the result from the integration of (1271) . 



3.2 A reformulation of the Kyner and Bennett approach 



In this section we provide the derivation of an alternative form of the equations of motion (|26[) . Starting from the kinematic 
representation (|19[) of the SPE we take the second derivative with respect to the time to obtain 



d 2 r 



di) ~dp rrtn + lf 



<ff df_ . 

at* at 



+ Rr r 



(28) 



In keeping with the Kyner & Bennett approach (equations d©J and (f8))- (|10|) 'l. every term is now reshaped into a form based 
on the true anomaly. 

For the scalar components r and r we need the first and second time derivatives of both the radius vector and the true 
anomaly. From (|2T|) we immediately find 



with n given by @, and subsequently 

where the perturbed gravitational parameter fi is given by 
" = n o (1 - if ffl3 - 

From the second of equations (|18l) follows that 



(29) 

(30) 
(31) 
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dr na 

^ = 7CT esm/ ' (32) 

which can be differentiated to yield 
d 2 r ecos/ 

Utt? ~ ^ ' (33) 
For the matrix terms we need the derivatives of the rotation matrix R with respect to true anomaly /. These are easily 
computed from (|22[) and with the help of (0 they can be written as 

dR dD , ^„dB , 

w = r ^ CB + (1 + ??)DC ^ (34) 

d 2 R 2d 2 D^„ „ .„ , dD„dB s2^^d 2 B ,„„, 

dp = ^ CB + 2r(1 + ^^ C ^ + (1 + ^ DC ^- (35) 

Direct use of these expressions clearly leads to the original K&B formulation. Our purpose here is instead to reformulate 

equations (|34[) and (|35[) through the introduction of the matrix operators Vf — d/df and V 2 — d 2 /df 2 . It is shown in 
Appendix [Althat 

V f = rV+(l + 77)H, (36) 

V 2 f = r 2 K + 2r(l + 77)N- (l+r?) 2 !, (37) 

where V, H, K and N are matrices defined in (]A5[) . (]A6[) . (|A7|) and IjAgp . respectively, and where I is the identity matrix of 
order three. Then 

— = © / R=[rV + (l + J] )H]R, (38) 

^ = V 2 R = [r 2 K + 2r(l + V )N-(l + r,) 2 ] R, (39) 
and 

(2?/R)r rin = [rV + (1 + 77) H] r, (40) 

(X>|R)r rtn = [r 2 K + 2r(l+77)N-(l + 77) 2 ]r. (41) 
Substituting equations pS]). ([30|>. ([32|>. (J33J), J40} and (gTJ) into ^28j we find 

f = n 2 (l - e 2 )^ [r 2 K + 2r (1 + 77) N - (1 + t?) 2 ] r + fi^^-r. (42) 

Notice that in performing the substitutions all traces of the first derivative DjR have been lost. Adding and subtracting the 
term (^/r 3 ) r, with fi from (|31|l . on the right hand side yields 



r + -^r = n 2 (l - e 2 ) ^ { [l - (1 + V ) 2 ] I + r 2 K + 2r (1 + 77) N} r. (43) 



A 4 2 
^3 r = 71 

This equation can now be put in the compact form 



.. . A 1 MPc 



(44) 



where p = a (1 — e j is the orbital semi-latus rectum, or parameter of the ellipse, and the matrix S has been defined as 
/ A B sin SI \ 

S = A B cos SI , (45) 
V C J 

with the constants 

A = 1-(1 + 7 ? ) 2 -t{t + 2(1 + 7?)cost:}, (46) 

B = 2r (1 + 77) sin i, (47) 

C = l-(l + 77) 2 . (48) 

In component form, with r = [x, y, z) T , the reformulated equations of motion of Kyner & Bennett read 
x + -^x = { [l - (1 + r/) 2 - t{t + 2 (1 + 77) cos i}] x + [2r (1 + 77) sin i sin SI] z} , 

V + ^y = { [l - (1 + ?7) 2 - r{r + 2 (1 + 77) cos i}] ?/+ [2r(l + 77) sin i cos SI] 2} , 



- 7T[1-(1 + ^) 2 ]- (49) 

Equations (|44|) . although highly simplified with respect to equations (J26J) — note the absence of matrix multiplications — 
still need to be supported by equation (|27p . since true anomaly appears implicitly through the node SI (see equation (©). 
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4 THE SPE IN THE MEAN ANOMALY 

This Section provides the formulation of the secu l arly p recessing ellipse with the time, or the mean anomaly, as the independent 
variable. We review the model due to lEscoball l|l965h and develop alternative formulations following essentially the same 
original line of development, but which are more concise for implementation. 

4.1 The Escobal formulation 

The procedure to obtain the equations of motion is well described in lEscobaH (| 19651 ). Here we give only the final result and the 
necessary auxiliary equation in case it is desired to fully integr ate the equatio ns of motion as a quick alternative to following 
the analytical update of the eccentric anomaly as envisioned in lEscoball (|l965h - The equations of motion are derived from the 
decomposition of the radius vector 

r = XP + YQ, (50) 

in the Laplace reference frame, where P is the Hermann- Jacobi-Laplace unit vector and Q is a unit vector normal to P in 
the direction of increasing anomaly. The unit vectors P and Q are themselves defined as 

P x = cos uj cos — sin w sin Q cos i, Q x — — sin uj cos fl — cos uj sin Q. cos i, 

P y — cos u sin Q + sin uj cos Q cos i, Q y = — sin uj sin Q + cos u) cos Q cos i, (51) 

P z = sin u; sin i, Q z = cos uj sin i. 

The acceleration is then 

f + -§r = 2XP + 2YQ + XP + YQ , (52) 
where the modified gravitational parameter jl is given by 

ft = n 2 a 3 , (53) 
and the first and second derivatives of P = (P x ,P y ,Pz) T and Q = (Qx,Q y ,Qz) T are given by 
P x = — tlPy + CjQ x , Q x = —ClQy — ljP x , 

Py = ClP X + djQy, Qy = O.Q X ~ dj Py , (54) 

P z =ujQ z , Qz=-diP z , 
and 

P x = - (O 2 + 6j 2 ) P x - 2Q.CuQ y , Q x = - (ti 2 + Cj 2 ) Q x + 2Q6jP v , 



P y = - (Q 2 +u 2 ) P y +2QujQ x , Q y = ~(n 2 +u 2 )Q y -2nujP x , (55) 

Pz = ~U 2 Pz, Qz = ~C0 2 Qz. 

Finally, X, X, Y and Y are defined either in terms of the eccentric anomaly E or in terms of the true anomaly / as 

na 2 na 
X = a (cos E — e) — r cos /, X = sm E — -^=^= sm /, 

VT ^ I 1 n ( 56 ) 
Y = aVT= ^sin E = r sin/, Y = ^- ^ V 2 cos E = na( ^ + C ° S/) , 

where the perturbed mean motion n is given by equation (|14|) . As in the true anomaly-based Kyner and Bennett formulation 
the number of first order differential equation to integrate is effectively seven, that is, six equations for the state vector and 
one equation (uncoupled with the state) for eccentric anomaly 

dE (57) 



dt 1 — e cos E ' 
or the true anomaly 
df n 



dt (i _ e 2)3/2 



(1 + ecos/) 2 . (58) 



Note that the equation for the anomaly is not redundant, but provides necessary information for the integration of the 
differential equations Q52p . 

4.2 A reformulation of the Escobal approach 

A simpler form of Escobal's equations of motion can be obtained by resolving Eq. (|52[) in inertial axes and substituting in it 
Eqs. (|54p . (|55p. (|51|l and (|56p for the true anomaly to obtain 

r+^r = -eA - {Br + A} cos / - {Cr + D} sin /, (59) 
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where the components of the vectors A, B, C and D are 
A x = 2na (CjP x + tlQy) , B x = (h 2 + Co 2 ) P x + 2(lCoQ y , 



vT 



e 



'6jP v - QQ X ) , B y = (ft 2 + Co 2 ) Py - 2flCoQ x 



A z = 2na uP z , B~ = Co 2 P, 



o2 



Inn ( 60 ) 
C x = -2DuP y + (fi 2 + Co 2 ) Q x , D x = [-VlPy + CjQ x ] , 

C y = 2Q,CoP x + (Cl 2 + Co 2 ) Q y , D y = 2na \ttP X + djQy] , 

Vl — e 

« .9/1 7-. 2nd . 

C z = u Q z , L>z = —==U)Q Z - 

vl — e 2 



Equations (|59p are written explicitly in the true anomaly, and this requires that for their numerical integration they be 
supplemented with equation (I58|l , 

Alternatively, one can choose to work with the eccentric anomaly, in which case the equations of motion read 

f+^r = eB'-|A'i + B'|cosS- j C" + D' i j sin E, (61) 

and the components of the several vector coefficient are 
e 2 )A. B' = aB, 



C = aVl-e 2 C, D' = aVl-e 2 D. 



(62) 



5 A NEW FORMULATION OF THE SPE IN THE MEAN ANOMALY M 

In this section we give the detailed procedure we have followed to obtain a new form of the equations of motion respresenting 
the SPE in the mean anomaly. Starting from the kinematic representation of SPE (|19[) and taking the second derivative with 
respect to the time we have 

d 2 r D d 2 r rtn dR dr rtn d 2 R 

■dF = ~dlc~2~ + ~dl~dT + ^ rrtn ' (63) 

which are to be used in connection with equations (|12|l , (|13[) , (|15|) and (|11|) that give the secular rates, the constant 7 and 
the evolution of the orbital elements with respect to the mean anomaly. 

The time derivatives of the rotation matrix R can be easily computed from (|22p recalling (|lip . As in the case of the K&B 
reformulation of Section ^. 21 the strategy is to define appropriate operators D t = d/dt and V 2 — d 2 /dt 2 acting on the rotation 
matrix R. In Appendix iBl it is shown that these operators can be expressed as 

V t = (w + /)H + nV, (64) 
V 2 = (Co + /) 2 H 2 +2tl(Co + /)N + 2 K + /H. (65) 

If we now recall that the use of the time as independent variable implies from the second of (|18p that 

e sin /, (66) 

(67) 



e cos / 



with fi given by (|53[) . we can write 

= TVtn, (68) 

(69) 



dv r tn 


nae sin / 


dt 


r\J\ — e 2 


d 2 r , 

' rtn 


e cos / 


dt 2 


= V r 3 r 



and substitute (|64[) . (|65l) . (16811 and (I69|l into (163[) . we get the equations of motion in the form 

r = fi^^r + 2 fmeSin/ \(Co + /)H + OVl r + \(Co + /) 2 H 2 + 2Vl(Co + f)N + ti 2 K + M r. (70) 
r J ry /i — g 2 J 

We now need expressions for the first and the second time derivatives of the true anomaly. These can be obtained from (1581) 
and put in the form 

2 



/ = rVl-e 2 ^, (71) 

/ = 2^L=/sin/. (72) 
ryl — e. 
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We thus easily see that in (|70[) the term with /H as a factor cancels with the term /H. If at the same time we add the modified 
Keplerian term fir/r 3 to both the left and right hand sides, the equations of motion become 

f + 4r = fi 1 + e ™ Sf r + 2 nae r ^L (CjH + ftV) r + [(* + /) 2 H 2 + 2Q(u + f)N + ft 2 K] r. (73) 

Considering that, on the basis of the relationship (|B17[) . the action of H 2 on r is simply to reverse its sign, 
H 2 r = H 2 Rr rtn = RKr rtn = -Rr rt „ = -r, (74) 
we can write 

,n2 u 2 -P _l + eCOs/ 

(/) H r = -fi^r = -n - 3 r, (75) 

and thus eliminate the first term on the right hand side of ()73|) with the term proportional to (/) 2 . We thus obtain 

r + ^r = ? n ° ^ sm / ^ H + ftV) r + [2ft(c0 + /)N + ft 2 K] r - (d; 2 + 2w/)r. (76) 

Collecting the terms depending on / and substituting / from (|71[) and eliminating n in favor of the modified area integral 



h = na 2 y / l~e 2 , (77) 
we obtain 

f+^r= (ft 2 K + 2ft^N-^ 2 l)r + 2^sin/(^H + ftV)r + 2-^ (ftN - wl) r. (78) 
If we now define 

E = ft 2 K + 2ft^N-^ 2 , (79) 

F = tOH + ftV, (80) 

G = ftN-tOI, (81) 

the equations of motion assume the more compact form 

r + ^-r = ( E + 2— sin/ F + ] r. (82) 

r j i p r r z J 

The matrices E and G are upper triangular, while matrix F is antisymetric. These three matrices are time-dependent through 
the longitude of the node, and are given explicitly as 

— (ft 2 + 2tluj cos i + lj 2 ^ 2ftcj sin i sin ft it) 
E(t) = | - (ft 2 + 2ftw cos i + Co 2 ) -2ftwsinicosft(i) ), (83) 

-Co 2 

— (ft + ujcosi) — Co sini cos ft (t) 

F(i) = | (tl + Cocosi) -wsinisinft(t) ] = -F T (<) , (84) 

w sin i cos ft (t) Co sin i sin ft (t) 

— (ft cos i + (jjj ft sin i sin ft (t) 
G(t) = | -(ftcosi + cO) -ft sin i cos ft (<) |. (85) 

-Co 

It appears odd that the longitude of the pericenter does not appear explicitly in the formulation. However, we need to 
recall that the explicit dependence on time is restricted to the position of the orbital plane, which depends on the inclination 
i and the node ft, only the latter of which is variable in our present model. 

The final and crucial step is to recognize from (|68|) that 

dr he . _ 

_ = — sin/ =r -r, (86) 
at p 



so that Eq. (|82|) can be reduced to the form 

(87) 



/i 

r + -j r 



_ 2,. „. „ 2ft „ 
E + - Or -r)F + — -G 

r r z 



The right-hand side of this equation is the explicit form of the disturbing acceleration s (r,r,t) appearing in equation {TJ. In 
component form, with r = (x,y,z) T , equation ()87[) reads 



x + -^x — — (ft 2 + 2QCo cos i + Co 2 ) x + (2C1C0 sin i sin ft) z + 2 XX 2 ZZ [— (ft + Co cos 1) y — (u sin i cos ft) 2] 
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uj sin i sin fl) z 



(88) 



(89) 



(90) 



This form of the equations of motion of the secularly precessing ellipse has the benefit of being a minimal and complete set 
of equations. In fact Eq. (|87[) has coefficients which are all explicit functions of the time and there is no longer any need for 
the supporting differential equation ((58")) for the anomaly /. 



6 DISCUSSION AND CONCLUSIONS 

The Kyner and Bennett and the Escobal equations of motion defining the intermediary orbit as an ellipse subject to secular 
motion of its angular elements have been reviewed. In both cases a more compact formulation has been developed and 
presented, which is better suited for implementation. Escobal's equations have also been reformulated and extended to allow 
the use of either the true or the eccentric anomaly, in addition to the mean anomaly. In fact, both these classical formulations 
are 'redundant, 'in that they include one more differential equation than the minimum of six associated with the three degrees 
of freedom of the problem. The reason for this extra equation is due to the particular choice of variables, which makes it 
necessary to propagate the perturbed anomaly, be it the eccentric or the true anomaly, along with the dynamical variables. 

The main contribution of the present work is the development of a novel formulation of the equations of motion of the 
secularly precessing ellipse that uses the time as the independent variable and requires no additional equation to account for 
the evolution of the anomaly. The final equation has the very compact form of a Two-Body equation of motion perturbed by 
a time-dependent acceleration containing three terms of degree and ±1 in the radius vector. The explicit dependence on 
the time is only due to the presence of the longitude of the node in the forcing terms. It should be noted that the supporting 
equation for the anomaly cannot be simplified away for true anomaly-based theories like Kyner and Bennett's. 

The time-wise approach developed here can be used when it is desired to numerically verify the analytical propagation of 
perturbations based on Kaula-type linear theories, where the nominal trajectory is a secularly precessing orbit. In particular, 
it can be used to verify the perturbation spectrum. It can also be used when analyzing first-order perturbation effects on 
orbital arcs, or ephemerides, estimated from observational data. In that case the secular rates of the angular elements to be 
used in the present, novel formulation can be estimated very accurately by numerically fitting the estimated orbital ephemeris. 
Clearly, since they are based on observational data, these rates include the effects of all acting secular perturbations. When 
applying the present formulation to evaluate the effects, for instance, of neglected sources of perturbations along the given 
orbit, it is therefore necessary to exclude all terms generating secular perturbations from the forcing acceleration g(r,r,t) 
appearing on the right-hand side of equation ([2]). In this particular application, if Encke's approach is adopted, no rectification 
is then needed, since the formulation is guaranteed not to generate any secular drift between the perturbed and the reference 
orbits. Applications of this novel formulation to the case of tidal perturbations will be the subject of a future contribution. 
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APPENDIX A: THE OPERATORS V F AND V% 

The operators 23/ and D 2 were defined in Section T3. 2 1 in terms of the matrices V, H, K and N. Comparison of equations (|34[) 
and (|35[1 with equations (|36p and (|37|l makes it clear that, since R is orthogonal, these matrices are defined as 

H = DC?R T , (A2) 
du 

d 2 D T 

K = ^CBR- (A3) 
Now since R T = B T C T D T , it is straightforward to show that 

(A5) 




— cos i — sin i cos Q 

•osi — sin i e 

sin i cos 17 sin i sin SI 



DC — B T C T D T =( cosi -sinisinfi I, (A6) 

du 




(A7) 

, n , R / — cos i sin i sin Q \ 

N = — C — B T C T D T =| -cosi -sinicosfi , (A8) 
d " du \ / 

where we have repeatedly used the fact that 

10 0). (A9) 

This equation, like (|A5[) . shows a familiar property of rotation matrices fcf. IParsI (|l98ll )). 

It remains to be shown that the action of the last operator DCB UU on the right hand side of equation (|35p on r rtn is 
equivalent to multiplication by — R. This follows immediately once it is realized that 

with 

/ \ 

(All) 





In fact, since Sr rtn = 0, we have 
d 2 B 

DC—— -r r tn = DCSr rtn - Rr rt „ = — Rr rtn = — r. (A12) 
du 2 

The full operator, analogous to V, H, etc., of course is defined as 

DC^-4r t = DCSR T I, (A13) 
du 

but for the present purposes it is expedient to use its restriction to vectors having a null third component, which justifies our 
definition of the operator in equation (|37[) . 



APPENDIX B: THE OPERATORS V T AND V% 

For the first derivatives of rotation matrix R (|22|) with respect to time we have 
dR dD„„ „„dB 

— = — CB + DC — . (Bl 

dt dt dt v ' 

The first term to the left hand side can be rewritten as 
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^CB = fi^D T DCB = fiVR, (B2) 
dt dQ K ' 



where V is the matrix already introduced in Eq. (|A5[I 
The second term of (IB It) can be reformulated as 



DC^ = «DC^ = (w + /)HR, (B3) 
dt du 



where H is given by Eq. (|A6[) . Collecting terms we have that the first derivatives of R is 

cm 

dt 



^ = [(a + /)H + nv]R. (B4) 



The second derivatives of R can be computed from previous equation that is 

^ = j t [(w + /)H + OV] R+ [(cj + /)H + fiV] 2 R. (B5) 
The first term of right hand side can be computed as 

j t [{u + /)H + nv] = /H + A(w + ,/)T, (B6) 
where we have substituted for 

H = ClT , (B7) 
with 

(0 sin i sin Q \ 

-sinicosfi . (B8) 

— sin i sin Q sin i cos $1 / 

Substituting and collecting terms in equation (|B5[) yields 

^ = [(w + /) 2 H 2 + 20(^ + /)N + Si 2 V 2 + / H] R, (B9) 

where we have used 

/ — 2cosi sinisinfi \ 

Q = HV + VH= -2cosi -sinicosn , (BIO) 

\ sin i sin — sin i cos Q / 

and 

Q + T = 2N, (Bit) 
with 

(— cos i sin i sin Q \ 

— cosi — sin i cos Q . (B12) 
/ 

It is now possible to define the operators Dt and T>\ as 

V t = (w + /)H + nV, (B13) 
O t 2 = (w + /) 2 H 2 +2Si(tO + /)N + 2 K + /H, (B14) 

where we have used the identity 

V 2 = K. (B15) 

Note that these operators are intended for application to a position vector only. A further relationship is useful in connection 
with the operator H 2 . If we insert the identity BB T after the factor C in equation (IA2|) and note that B T (dB/du) = V, it 
follows that 

H = RVR T . (B16) 
Then, taking the square of both sides and multiplying by R on the right we find 

RK. (B17) 



This paper has been typeset from a T^X/ ETgX file prepared by the author. 
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